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CNSFV Code Development, 

Virtual Zone Navier-Stokes Computations of 
Oscillating Control Surfaces and 
Computational Support of the Laminar Flow Supersonic 

Wind Tunnel 

Abstract 

The work performed during the past year on this cooperative agreement covered 
two major areas and two lesser ones. The two major items included further develop- 
ment and validation of the CNSFV code and providing computational support for 
the Laminar Flow Supersonic Wind Tunnel (LFSWT). The two lesser items involve 
a Navier-Stokes simulation of an oscillating control surface at transonic speeds and 
improving the basic algorithm used in the CNSFV code for faster convergence rates 
and more robustness. The work done in all four areas is in support of the High 
Speed Research Program at NASA Ames Research Center. 

Introduction 

The numerical simulation of the Navier-Stokes equations for complex configua- 
tions at realistic flight conditions is still limited by several problems. Included are 
the lack of adequate grid resolution, robust and efficient flow solvers, transition pre- 
diction techniques, and turbulence models. The work covered in this report involves 
the first three items only. The grid resolution and robust flow solver problems are 
included in the CNSFV code development work and the transition prediction tech- 
nique is included in the work involving flow analysis of the laminar flow supersonic 
wind tunnel currently under development at NASA Ames. The flow analysis is by 
numerical simulation of the Navier-Stokes equations with the CNSFV code. In the 
following sections the work performed during the past year in the four different 
categories are described in more detail. 

CNSFV Code Development 

The CNSFV code is a cell-centered finite volume version of the finite differ- 
ence code, CNS. The reason for developing the finite volume code was to facilitate 
the implementation of conservative zonal interface boundary condition. The finite 
difference form is not practical for conservative interfacing. The CNSFV code was 
originally developed three years ago (see ref. 1) and has been under continuous 
improvement since its inception. During the past year further improvements have 
been made. These include simplifying the type and amount of input data required, 
implementing general boundary conditions, increasing the code efficiency in terms 
of vectorization and algorithmic improvements, and developing special zoning capa- 
bilities (called virtual zones) to ease the problem of generating grids with multiple 
zones about complex and dynamic aerodynamic configurations. An example of an 
application of the virtual zone technology is given in reference 2, where a complete 
wing-body configuration with control surfaces is simulated with the CNSFV code. 

The diagonal ADI algorithm in the original finite difference code is a fast and 
robust scheme. However when the code was converted to finite volume form a sharp 
drop in the maximum allowable stable time step was noticed. Various types of local 
time step scalings were tried and the allowable times steps improved dramatically 



but the time steps originally possible with the finite difference form were never 
achieved. The scaled allowable time steps were still too small however and the 
diagonal ADI scheme was replaced with the lower-upper symmetric Gauss Seidel 
(LU-SGS) scheme. This scheme is unconditional stable and arbitrarily large time 
steps can now be used. Most of the the work on implementing and validating the 
LU-SGS scheme in the CNSFV code is reported in an AIAA paper presented at the 
AIAA Fluid Dynamics Conference, ref.3. This paper is included as Appendix A to 
which the reader is referred for a complete discussion. 

A preliminary version of an user’s manual (ref. 4) for the CNSFV code has 
been written and is presently undergoing user testing. This work is still underway 
and expected to be completed by the end of this year. 

Computational Support for the LFSWT 

For the past 18 months computational support has been provided for the Lam- 
inar Flow Supersonic Wind Tunnel (LFSWT). The objective of the effort is to 
develop computational tools so that the design of a test model and its placement 
within the test section of the LFSWT can be verified by numerical simulation of 
the Navier-Stokes equations before the model is constructed. For transition stud- 
ies in the supersonic Mach regime it is important to know the extent of clean and 
undisturbed flow over the test model. 

For the simulation, modified versions of the Upwind Parabolized Navier-Stokes 
(UPS, ref. 5) and the Compressible Navier-Stokes, Finite Volume (CNSFV) codes 
were used to solve the thin-layer Navier-Stokes equations for the laminar flow about 
the test model inside the LFSWT. The surface and flow field grids were generated 
with GRIDGEN (ref. 6). The faster UPS code was used for the higher Mach 
numbers investigated and the multi-zonal CNSFV code for the lower supersonic 
Mach regime. 

Computations have been performed of flow fields about a NACA64A010 wing 
with a 70° leading edge sweep mounted on the top wall of the LFSWT for inviscid 
and viscous flows with the UPS and CNSFV codes, respectively. Various other 
model locations were studied to verify that the top wall mounted position provides 
for the largest extent of undisturbed flow on the model. Figure 1 shows the inviscid 
shock pattern obtained with the UPS code. Shown are the impinging shocks on the 
tunnel walls as well as on the model itself. The undisturbed region on the model is 
the triangular region in front of the reflected shock wave impinging on the model. 
The flowfield behind the impinging shock wave is no longer undisturbed and, hence, 
that part of the model behind the shock is useless for any natural transition study. 

Another series of computations were carried out on the NACA64A010 wing 
with four different leading edge sweep angles. The mean flow results were used to 
validate a Parabolized Stability Equation (PSE) code being developed by another 
group. In addition a full scale sized portion of the F16XL2 passive glove is presently 
being gridded and a flow solution is being obtained in preparation for the wind 
tunnel tests scheduled for early 1994. The results of the computations axe being 
used to design the wind tunnel model of the F16XL2 wing. The wind tunnel tests 
are designed to validate the wind tunnel with flight tests. 



An accurate prediction of the flow field inside the LFSW with various test 
models is important for maximizing the usefulness of the tunnel, especially when 
relatively small test sections are considered. Numerical simulations can also be used 
for designing tunnel modification and innovative passive and active tunnel devices 
to minimize the impact of reflected waves on the test model. This is a cost-effective 
means of increasing the usable size of the LFSWT. 

Virtual Zone Navier-Stokes Computations of Oscillating Control Surfaces 

Another area of effort conducted during the past year was to implement the 
virtual zone concept into a time accurate finite difference code, ENSAERO (ref. 7), 
for application to an oscillating control surface mounted on a clipped delta wing at 
transonic speeds. For a time accurate computation it is essential that the search 
procedure to find the interpolation coefficients for the inter-zonal communication be 
at least as efficient as the flow solver. Otherwise the code is not practical enough for 
routine use. Much of the work involved developing more efficient search procedures. 
The results of this effort were presented at the AIAA 11th CFD Meeting (ref. 8). 
The paper is included in this report as Appendix B. 

Algorithm Development 

As mentioned the CNSFV code development section, the LU-SGS scheme is un- 
conditionally stable. However it does suffer slower convergence rates with increasing 
Reynolds number for viscous dominated flows. For this reason several modifications 
of the basic LU-SGS scheme were implemented. One modification succeeded in im- 
proving the convergence rates by as much as a factor of four. A typical example of 
the improved convergence rate for a compressible flat plate boundary layer flow at 
Mach 2 is shown in Figure 2. This work is still underway and the complete results 
will be reported at a future date (AIAA 25th Fluid Dynamics Conference in July 
1994). 
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Figure 1. Shock pattern on a NACA64A010 wing in the LFSWT. 
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Figure 2. Comparison of convergence rates between standard and 
modified LU-SGS schemes for a flat plate boundary flow. 
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Abstract 


The LU-SGS (lower upper symmetric Gauss Seidel) 
algorithm has been implemented into the Compressible 
Navier- Stokes, Finite Volume (CNSFV) code and validated 
with a multizonal Navier-Stokes simulation of a transonic 
turbulent flow around an Onera M6 transport wing. The 
convergence rate and robustness of the code have been im- 
proved and the computational cost has been reduced by 
at least a factor of 2 over the diagonal Beam- Warming 
scheme. 


Introduction 


The numerical simulation of the Navier-Stokes equa- 
tions about complex and realistic aerodynamic configu- 
rations with a structured grid requires the use of zonal 
methods. A popular and fairly efficient numerical scheme 
used to solve the Navier-Stokes equations is the diagonal 
implicit Beam- Warming algorithm [1,2]. A finite volume 
multi- zonal Navier-Stokes code, CNSFV, was developed 
using this scheme. However the Beam- Warming scheme 
uses approximate factorization to simplify the matrix in- 
version process. As a result the convergence properties of 
the scheme are dependent on the time step and time step 
scaling chosen, and much user input is required to deter- 
mine the optimum time step and time step scaling. For a 
numerical simulation of a realistic aerodynamic configura- 
tion, several dozen zones may be required. Choosing the 
optimum time step and scaling for each of these zones be- 
comes a tedious and time consuming process. For reasons 
not understood yet, the finite volume formulation of the 
diagonal Beam- Warming scheme suffers a severe time step 
restriction. Typically for constant time steps, the largest 
CFL number that could be exercised were less than ten 
and much less than that for the initial transients. With a 
judicious time step scaling, the maximum CFL can be as 
high as 50 to 75. The time step restriction is not so severe 
with the finite difference formulation of the diagonal Beam- 
Warming scheme, where the maximum CFL numbers can 
be as high as 500. The time step restriction for the finite 
volume form of the diagonal Beam- Warming scheme is se- 
vere enough that the time step is too small to be practical 
for an unsteady Navier-Stokes code. 


* Senior Research Scientist, MCAT Institute, San 
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A numerical scheme without the above mentioned 
time step restrictions is the LU-SGS (lower upper symmet- 
ric Gauss Seidel) algorithm [3,4]. This scheme has been im- 
plemented into the CNSFV code and validated with a mul- 
tizonal Navier-Stokes simulation of a transonic flow around 
an Onera M6 transport wing. 


Navier-Stokes Equations 

The three-dimensional thin-layer Navier-Stokes equa- 
tions in strong conservation law form in curvilinear coor- 
dinates are 

drVQ + dtF + dnG- 
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The contravariant velocity components are defined as 


U = V((t+ + £ yV + (;W), 
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and the viscous fluxes are given by 
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The finite volume metrics represent the cell face area 
normals in each of the curvilinear coordinates (£,r/,£). 
They are related to the metrics introduced in equations 
(1 - 5) as follows 


The pressure is given by the equation of state 


p = (7 - l)(e - p ( u 2 + v 2 + tv 2 )/2) (6) 

where 7 is the ratio of specific heats. The sound speed 
is denoted by a. The nondimensional parameters are the 
Reynolds number Re and the Prandtl number Pr. The 
coefficient of viscosity fi and thermal conductivity k axe 
decomposed into laminar and turbulent contributions as 
follows: 


^ = m + ih 
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where Pri and Pr t are the laminar and turbulent Prandtl 
numbers and k = \x with the nondimensionalizationused in 
this paper. The standard Baldwin-Lomax turbulent eddy 
viscosity model [5] is chosen for this study. 

The nondimensional parameters chosen for this code 
are the same as those in the CNS code [6]. Denoting di- 
mensional quantities with a (), the normalizing parameters 
are the frees tream density poo, the freestream sound speed 
doo, the freestream viscosity coefficient ft oo, and a charac- 
teristic length l. 

The metrics used above have a different meaning for a 
finite volume formulation compared to the finite difference 
formulation of [6]. Referring to a typical finite volume cell 
as shown in figure 1, the finite volume metrics are defined 
(see, for example, Vinokur [7]) as 
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where x* = x,y,z for i = 1,2,3. The volume of the com- 
putational cell is given by 


V = ^[( r 4 - r 2 ) x (r 3 - ri) * (r 7 - r 3 ) 
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and is the finite volume equivalent of the inverse Jacobian 
of the coordinate transformation in the finite difference 
formulation of [6]. 

Numerical Method 

The governing equations are integrated in time for 
both steady and time accurate calculations. The unfac- 
tored linear implicit scheme is obtained by linearizing the 
flux vectors about the previous time and dropping second 

and higher order terms. The resulting scheme in finite 

volume form is given by 
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where the residual R n is 
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The convective flux jacobians A, B, C and the viscous flux 
jacobians L, M and N are defined in the appendix of [1], 

Solving the above equation set by a direct matrix in- 
version is still not practical for three dimensional problems. 
However there are several indirect or approximate methods 
available, including the diagonal Beam- Warming scheme 
and the lower- upper (LU) factorized scheme of Yoon and 
Jameson. In a previous study [8,9], the diagonal Beam- 
Warming scheme was the basic flow solver algorithm for 
the multi-zonal compressible Navier- Stokes finite volume 
code, CNSFV. While the scheme worked well for finite dif- 
ference codes, the performance deteriorated for the finite 
volume code. The cause of the deterioration is not known, 
but it is not isolated to the CNSFV code; several other 
finite volume codes using ADI schemes seem to suffer from 
the same type of deterioration. Because of this problem, a 
LU scheme which does not seem to suffer the same problem 
as the ADI schemes is incorporated into the CNSFV code. 
A brief description of both the diagonal Beam- Warming 
and the LU scheme is presented. 


Each of the factors of the implicit operator of equa- 
tion (12) has an artificial dissipation term added to sta- 
bilize the central difference operator. The added term is 
based on Jameson’s nonlinear second and fourth order dis- 
sipation and for the ^-operator takes the form 
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Diagonal Beam- Warming Algorithm 

With the use of approximate factorization and diag- 
onalization of the flux jacobian matrices, a scalar pentadi- 
agonal algorithm [2] can be derived from eqn (10) as 

Trfl-h V“ 1 A^A^]iV[/+ V- 1 h6 r ,A„)P[I + V^liS^} 

■T < -'AQ n = R n (12) 

where 8 $ is a central difference operator and A Q n = 

Q n+1 — Q n with Q n+1 = Q(t n + h). The viscous terms 
are not included in the implicit side. The artificial dissi- 
pation is included in both sides and is derived below. 

The inviscid flux jacobians are diagonalized as fol- 
lows: 

OqF = A = T^Tf 1 
8qG = B = TjjArjT" 1 
OqH = C = T c A < T~ l 

The , Tt,, T( are the eigenvector matrices of A, B, C, 
respectively with A^, A^A^ as the respective eigenvalues. 

We also have 


where k 2 , «4 are constants of o(l), and A^, V* are the for- 
ward and backward difference operators. &j+k is a modi- 
fied spectral radius defined as 
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and where cell centered surface areas are used, e.g. 


The modified form of the spectral radius, equation 
(17), is suggested by Turkel [10] to account for large aspect 
ratio computational cells as for example in a viscous layer. 
Similar dissipation terms are obtained for the rj— and £- 
operators. The dissipation terms added to the right hand 
side of equation (10) are identical to those given above 
except that A Q n is replaced by Q n . 


LU-SGS Scheme 
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The unfactored scheme of eqn (10) is given in terms 
of central differences for the implicit operator. It can also 
be represented in terms of upwinded differences. Dropping 



the viscous terms in the implicit operator, equation (10) is 
in terms of upwind differences as follows: 


{I+V- 1 h(6~A++6+A-+6;B++6+B~+6;C++6+C-)] 

•A Q n = -V- 1 hR n (20) 

where 8^ and <5^* are the forward and backward difference 

operators. Similarly, A + and A~ are the Jacobian matrices 
which contain non-negative and non-positive eigenvalues, 
respectively. 

The Yoon and Jameson version of the LU scheme can 
be obtained from the above equation by a simple reordering 
of the matrix elements and approximately factoring into 

two matrices. Define D ) L, and U to be matrices which 
contain the diagonal, sub-diagonal and super-diagonal el- 
ements of the implicit operator of equation (20), respec- 
tively. 

[D + L 4- U]AQ n = —V~ l hR n 
which can be factored into 


[D + L\D - 1 [D + U]AQ n = - y- 3 hR n 

Redefine L = D + L and U = D + U to yield the final form 
of the LU-SGS scheme. 


LD~ l UAQ n = -V " 1 hR n (21) 

where 

L = I 4- V- l h(6£ A + + + 6~C+ - A~ - B~ - C") 

D = I + y _1 /i(A+ - A- + B+ - B~ + C + - C") (22) 

U = I 4- V-' h(6+A- 4 8+B- + 6+C- + A + 4- 5+ + C + ) 

A variety of LU-SGS schemes can be obtained by differ- 
ent choices of the numerical dissipation models and Jaco- 
bian matrices of the flux vectors. In this study we use the 
same artificial dissipation described for the diagonal Beam- 
Warming scheme. For robustness and to ensure that the 
scheme converges to a steady state, the matrices should 
be diagonally dominant. To ensure diagonal dominance, 
the Jacobian matrices can be constructed so that 4- ma- 
trices have nonnegative eigenvalues and — matrices have 
nonpositive eigenvalues. For example, the diagonalization 
used for the Beam- Warming scheme can be used to obtain 
the zb matrices. 


A* = T^Af Tf' 


B* = T.A^T- 1 (23) 

C* = TcAfT ^ 1 

Another method to obtain diagonal dominance is to con- 
struct approximate Jacobian matrices: 


A ± = [A±p(A)I }/ 2 

B ± = \B±p{B)I}l2 (24) 

C ± = [C±p{C)I }/ 2 

where /5(A) = maa;[|A(A)|] and represent a spectral 

radius of the Jacobian matrix A with the eigenvalues 
A(A). These eigenvalues are, e.g., A(A) = [£/, 17, 1/, U ± 

a yJ&i 4- sj + s\\j, where the metric terms are again de- 
fined at cell centers by Eq. (19). Other methods of increas- 
ing the diagonal dominance of the LU operator include the 
addition some approximation of the viscous terms and the 
artificial dissipation in the implicit operator [11]. 

The inversion of equation (21) is done in three steps. 
The block inversion along the diagonal is eliminated if the 
approximate Jacobian of Eqn (24) are used instead of Eqn 
(23). A Newton-type of iteration is obtained simply by set- 
ting h = oo. Eliminating the time step from the algorithm 
provides the practical advantage of side-stepping the need 
to find an optimal CFL number or the time step scaling to 
reduce the overall computational effort to achieve steady 
state solutions. As mentioned previously, this is an im- 
portant consideration for a multizonal code. If first-order 
one-sided differences are used, Eq. (22) reduces to 


L - pi - ( - Bj k _ i , - j 


D - pi (25) 

U = Pi + lfJM 4- #>,*+!,/ + £,-,*,/+ 1 

where 

p = /5(A) 4- p(B) 4- p(C) 

This algorithm requires only scalar diagonal inversions 
since the diagonal of L or U = D. The true Jacobians 
matrices of Eq. (23) may permit better convergence rates, 
but require block diagonal inversions with approximately 
twice the computational effort per iteration. This study 
considers the scalar version without the viscous or artifi- 
cial dissipation only. 

The scheme is completely vectorizable on j -f k 4- / = 
constant oblique planes of sweep, Fig. 2. This can be 
achieved by reordering the three dimensional arrays into 
two-dimensional arrays as follows: 


4 



P(i P ,p) = P U,k,l) 

where i p is the serial number of the oblique plane of sweep 
and p is the address of point (j, fc,/) in that plane. The 
specific value of p in terms of (j, fc, /) is 


p — j + k + l — 2 

and the number of points in the plane p is 
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The maximum number of points in the planes is n Pimax = 
n p (pi) } where pi — (p ma x + l)/2. The vector lengths 

change from 1 at the beginning of the inversion sweep to a 
maximum of n p ^ max at the center of the sweep and back to 
1 at the end. While n p ^ max can be large, the vector lengths 
at the beginning and end of the sweep are very small and 
inhibit efficient vectorization. Nevertheless for reasonably- 
sized zones, n Ptmax is sufficiently large so that good vec- 
torization can be achieved. Other reordering sequences are 
possible, but additional sweeps will then be necessary and 
the overall performance of the LU-SGS inversion is less. In 
addition, the oblique planes of sweep allow for very efficient 
autotasking procedures on multi-processor computers; the 
efficiency is less with other sweep directions. 


Boundary Conditions and Zonal Interfaces 

To complete the equation set boundary conditions 
must be specified. With the use of curvilinear coordi- 
nates the physical boundaries have been mapped into com- 
putational boundaries, which simplifies the application of 
boundary conditions. The boundary conditions to be im- 
plemented for external viscous or inviscid flows include (1) 
inflow or far field, (2) outflow, (.3) inviscid and (4) viscous 


impermeable wall, and (5) symmetry conditions. For ex- 
ternal three-dimensional flow fields about closed bodies, 
the topology of the grid usually introduces (6) grid singu- 
larities which require special boundary conditions. The use 
of zonal methods can avoid the generation of grid singu- 
larities, but requires (7) special zonal interface boundary 
conditions. For compressible flows these zonal boundary 
conditions must be conservative to maintain global con- 
servation. In this study we consider nonconservative in- 
terfaces. However, conservative interzonal communication 
has been developed and reported in a previous study [8]. 

A patched zonal procedure is used for the grid sys- 
tem. The zonal interfaces are general three-dimensional 
surfaces. The zoning procedure has a general capability 
so that a face of one zone can be in contact with several 
other zones. Information is transfered between zones with 
a bilinear interpolation procedure. The search method for 
determining the interpolation coefficients is automatic in 
that no user input is required other than identifying the 
zonal faces in contact with each other. The grid topology 
of the various zones is arbitrary and the zonal interfaces 
are not, in general, mesh continuous. In this study we use 
the ghost cell concept to obtain boundary and interface 
fluxes. Although there is no grid overlap at the zonal in- 
terfaces, the cell- centered conserved variables are reflected 
across surface and zonal boundaries and are determined by 
interpolating from the the adjacent zonal values. To deter- 
mine the dissipative fluxes for the fourth-order dissipation 
requires two ghost cells. However, here we simply reduce 
the fourth-order dissipation to second-order so that only 
one ghost cell is needed. 

Results 

The test case chosen to validate the LU-SGS scheme 
and the code is the attached turbulent transonic flow over 
an ONERA M6 wing. This case has good experimen- 
tal data [12] available and has been used extensively for 
Navier-Stokes code validation. The reason for choosing 
an attached flow case instead of a separated flow is that 
the emphasis of this study is to validate the flow solver and 
not the turbulence model. For attached flows the Baldwin- 
Lorn ax model [5] is sufficiently accurate. The flow condi- 
tons are = 0.84, a = 3.06°, Re mac = 11.72 x 10 6 . The 
wing surfaces are adiabatic and the fluid is a perfect gas 
with 7 = 1.4. The computational grid consists of a C-0 
mesh split into four zones. There are two zones around the 
wing, one each on the upper and lower surface of the wing. 
Similarly, the wake region is split into an upper and lower 
grids. This grid, in the single zone version, seems to be the 
standard grid used for various validations of Navier-Stokes 
codes, see, e.g., ref. 13. Figure 3 shows a partial view of 
the grids and the zones. The total grid size is 193 x 34 x 49 
points, and the individual zones consist of 73 x 34 x 49 and 
25 x 34 x 49 points for the wing and wake zones, respec- 
tively. The wall normal grid spacing at the surface is on 
the order of t/ + = 4. This is sufficient to properly resolve 
the boundary layer. 

The converged results in terms of pressure contours 
on the upper surface of the wing are compared with ex- 
perimental data [12] are shown in Fig. 4. The agreement 
between computation and experiment is good for the four 
span stations shown. These results are comparable to the 
results reported in ref. 13. The numerical predicted shocks 
are slightly smeared as compared to experiment. This is to 
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be expected since the numerical scheme uses central differ- 
encing with nonlinear scalar dissipation. Increased mesh 
resolution or an upwind scheme will sharpen the shock 
wave structure. 

Figure 5 presents the performance comparison be- 
tween the diagonal Beam- Warming scheme and the LU- 
SGS scheme. These views show the convergence rate for 
the upper wing zone with the two schemes. The first view 
shows that the LU-SGS scheme is not that much faster 
than the diagonal Beam- Warming scheme in terms of iter- 
ations required to reach convergence. However, in terms of 
cpu times, the LU-SGS scheme is at least twice as efficient, 
as shown in the second view. 

Most of the gains in efficiency with the LU-SGS 
scheme is in the reduced operation count. Table 1 shows 
the unit cpu time ( f.tsec / gridpoint/ iteration) for both 
schemes, as measured on the NAS A- Ames Cray C90 com- 
puter. The unit cpu time for the right hand side (RHS) in- 
cludes the calculation required for the RHS fluxes, bound- 
ary conditions, interface interpolation, swapping data 
in/out of machine core for each zone at each time step, 
writeout of convergence data, and other miscellaneous op- 
erations. This unit time is the same for both schemes at 
4.5^sec. The left hand side (implicit) unit cpu time is sub- 
stantially different. The times are 5.9 fisec for the diagonal 
Beam- Warming and 2.2 fisec for the LU-SGS scheme. The 
total times are 6.7 fisec and 10.4^sec for the LU-SGS and 
diagonal Beam- Warming schemes, respectively. Further- 
more, the diagonal Beam- Warming scheme required user 
intervention to obtain the optimum time step size as well 
as the time step scaling. The LU-SGS scheme uses no time 
step scaling and an infinitely large time step as explained 
in LU-SGS section above. 

Conclusions 

A multi-zonal compressible Navier-Stokes code has 
been improved by replacing the implicit diagonal Beam- 
Warming algorithm with the lower-upper symmetric Gauss 
Seidel scheme. With the new scheme the code is now much 
more robust, requires no user intervention to determine the 
optimum time step or time step scaling, converges faster, 
and requires only half the cpu time to obtain the same level 
of convergence. Most of the gain in efficiency is due to the 
lower operation count required by the LU-SGS algorithm. 
The algorithm and code have been validated with a tur- 
bulent simulation of an attached transonic flow around an 
Onera M6 wing. 
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Scheme 

CNSFV - diag BW 


Total 

10.4m sec (340 mflop) 

6.7psec (330 mflop)* 

LHS 

5.9fisec 

2.2psec 


xiinv 

lusgs 


etainv 

(oblique 


zetainv 

planes of 


vpenta, etc 

sweeps) 

RHS, 

4.5psec 

4.5psec 

etc. 




* average vector length ~70, optimum length ~128 


Table 1. Perftrace performance comparison of ding- 
onal Beam- Warming and LU-SGS schemes on Cray C90 
computer, (unit cpu time /tsec/grid point/iteration) 
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Fig. 1. Finite volume cell nomenclature. 


Fig. 2. Oblique plane of sweep for vectorization. 



Fig. 3. Four zone C-0 Grid for the ON ERA MG 
wing; only the two zones on the upper surface of wing and 
wake are shown. 
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CPU TIME (RELATIVE UNITS) 


Fig. 5. Comparison of the convergence rates ver- 
sus iteration and cpu times between the diagonal Beam- 
Warming and LU-SGS schemes. 
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Abstract 

A new zoning method called “virtual zones” has been 
developed for application to an unsteady finite difference 
Navier- Stokes code. The virtual zoning method simplifies 
the zoning and gridding of complex configurations for use 
with patched multi-zone flow codes. An existing interpo- 
lation method has been extensively modified to bring the 
run time for the interpolation procedure down to the same 
level as for the flow solver. Unsteady Navier-Stokes com- 
putations have been performed for transonic flow over a 
clipped delta wing with an oscillating control surface. The 
computed unsteady pressure and response characteristics 
of the control-surface motion compare well with experi- 
mental data. 

Introduction • 

Present transport aircraft as well as highly maneu- 
verable fighter aircraft are often subject to unsteady aero- 
dynamics. In this unsteady environment aircraft designers 
utilize active controls to achieve controllability and safety 
of the aircraft. Active control can also be used to suppress 
transonic flutter characteristics of high aspect ratio wings 
and thus reduce the structural weight to achieve more ef- 
ficient flight conditions. 

In the transonic flow regime active controls have a 
pronounced effect on the aerodynamic and aeroelastic per- 
formance of a wing. This effect can be used to improve the 
airplane performance by improved design of the active con- 
trol surfaces. To do this successfully requires accurately 
predicting the aerodynamic and aeroelastic performance 
of a wing. Experimental prediction of unsteady aerody- 
namic and aeroelastic performance is costly, risky and time 
consuming; numerical simulation of the unsteady Navier- 
Stokes equations is a much more cost effective alternative 
for predicting the performance of an active control surface. 

The physics of unsteady transonic flow around a con- 
trol surface has been simulated with small disturbance the- 
ory [1,2]. Unsteady Navier-Stokes simulations have been 
performed in two dimensions [3,4]. A more recent study [5] 
conducted a three dimensional simulation of the unsteady 
thin-layer Navier-Stokes equations on the flow field sur- 
rounding a wing with a forced oscillating control surface. 
In that study an unsteady Navier-Stokes code, ENSAERO, 
was extended to simulate unsteady flows over a rigid wing 
with an oscillating trailing-edge flap. 
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The numerical simulation of the unsteady Navier- 
Stokes equations about complex and realistic aerodynamic 
configurations requires the use of zonal methods. In this 
method the overall flow field domain is subdivided into 
smaller blocks or zones. In each of these zones the flow field 
is solved independently of the other zones. The boundary 
data for each zone is provided by the neighboring zones. A 
major difficulty of the zonal methods applied to oscillating 
control surfaces has been how to account for the variable 
exposure of the ends of the control surfaces to the flow 
field. 

In Ref. 5 an algebraic grid generation technique was 
incorporated into the ENSAERO code. The grid moved 
at every time step to follow the deflection of the flap. The 
small unsteady deflections were handled using a sheared 
single mesh. The large mean deflections were handled us- 
ing a zonal method [6]. The use of the single sheared grid 
did not permit the exact simulation of the unsteady flap- 
wing geometry. A gap had to be introduced between the 
ends of the flap and wing to allow sufficient space for the 
moving sheared mesh. The gap compromised the geom- 
etry and numerical predictions of the oscillating control 
surface flow field. The purpose of the present study is to 
rectify that compromise through the use of a new zoning 
technique called “virtual” zones. 

The virtual zoning method, first implemented in a 
multizone finite volume code, CNSFV, [7], has been mod- 
ified for application to the unsteady finite difference code, 
ENSAERO. The main purpose of these zones is to convert, 
for example, a solid wall boundary condition into an inter- 
face condition. The interface conditions are required for 
the interzonal communication. In a multi-zonal code- the 
virtual zones are treated like real zones as far as boundary 
and interface conditions are concerned, however, no flow 
field computations are done within these zones. Hence, 
the name “virtual” zone is appropriate. 

In addition to the introduction of the virtual zones, 
it is necessary to speed up the process of determining the 
interpolation coefficients required for the interzonal com- 
munication if the unsteady Navier-Stokes simulation is to 
be practical. 

The present study considers the transonic vortical 
flow over a clipped delta wing. A view of the wing and the 
control surface is shown in Fig. 1. Unsteady Navier-Stokes 
computations for the clean wing were reported in Ref. 8. 
The forced oscillating control surface computations with 
the single zone sheared mesh were presented in Ref. 5. 

Numerical Method 

A brief description of the governing equations, grid 
and zonal system, and boundary conditions is given in this 
section. Most of the attention will be focussed on the vir- 
tual zoning concept and the unsteady interpolation pro- 
cedure since these two items were crucial to the successful 
and practical application of the Navier-Stokes equations to 
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oscillating control surfaces. 


Governing Equations and Discretization 

The governing equations are the Reynolds- averaged 
thin-layer Navier-Stokes equations. The laminar viscos- 
ity is taken from the freestream laminar viscosity and is 
assumed to be constant for the transonic flow considered 
in this study. The turbulent viscosity is obtained with the 
Baldwin- Lomax algebraic eddy viscosity model [9] with the 
Degani-SchifF modification [10] to properly handle the lead- 
ing edge separation as well as the control flap vortical flow. 

The numerical algorithm, time dependent metrics of 
the curvilinear coordinate system, and performance char- 
acteristics of the ENSAERO code have been described 
previously and will not be repeated here. The interested 
reader is referred to Ref. 5. 


Control Surface Grid and Zones 

The primary focus of the present study is to demon- 
strate the feasibility of using dynamic zones for the os- 
cillating control surface case rather than minimizing the 
cpu run time for the case to be presented. Hence the flow 
domain was split into only three real zones. Each of the 
zones consists of a C-H topology. The two zonal bound- 
aries were placed at the span stations located at the ends 
of the control surface. Four additional virtual zones were 
placed in the two cuts separating the flap and the wing. 
One pair of the virtual zones remains fixed with the wing 
and the other pair is fixed with the flap and moves with 
the flap during the control surface motion. Figure 2 shows 
the seven zones used in this study. 

The C-H grid around a deflected control surface can 
be obtained in two ways. One is to shear every grid line 
normal to the control surface with the local deflection. The 
other is to algebraically regenerate the entire C-H grid with 
the control surface deflected at every time step. A previous 
study [5], showed that the computed surface pressures did 
not show any significant differences between the two meth- 
ods. Therefore, for this study, the grids around the control 
flap were regenerated with the simpler shearing method. 


Virtual Zone 

The zoning capability of the CNSFV code [6] allows 
the possibility of a single face of a zone to interact with sev- 
eral other zones. The procedure of determining the inter- 
polation coefficients is automatic in that no additional in- 
formation is required other than identifying the faces that 
are in contact with each other. To further extend the flex- 
ibility of the zoning method for the case of control surface 
aerodynamics, the idea of “virtual” zones was introduced 
in Ref. 7. 

Virtual zones are zones of zero thickness (for a fi- 
nite volume formulation) which serve to transfer solid wall 
(or other) boundary conditions to an interface condition. 
Thus multiple boundary conditions can be imposed on a 
block face with the same flexibility as an interface condi- 
tion. Virtual zones also decouple the process of volume 
grid zoning from the surface grid patches which define 
the aerodynamic configuration under study. Surface grid 
patches are required in order to impose the proper bound- 


ary conditions. The volume grid zones should be set up 
to obtain the proper mesh qualities required for numeri- 
cal accuracy. Another advantage of the decoupling is that 
much fewer zones are now needed, thus easing the effort 
and time required to generate the grids about complex and 
realistic configurations. 

Perhaps the simplest way to explain the virtual zon- 
ing concept is through the use of a generic wing/flap con- 
figuration represented by two cylinders sliding past each 
other as shown in Fig. 3. In Fig. 3a(i) the two cylinders 
are in contact and aligned. A convential zoning scheme 
has no difficulties for this configuration. In Fig. 3a(iii) the 
two cylinders are completely separated from each other and 
again a convential zoning scheme is adequate. However a 
convential zoning scheme may not be possible at the in- 
stance where the two cylinders are in partial contact as in 
Fig. 3a(ii). The topology of the grids and zones in these 
three instances is not the same, and thus it is not possi- 
ble to conveniently use the convential zoning procedure for 
such dynamic configurations. 

On the other hand, as Fig. 3b illustrates, with the 
virtual zone concept the topology of the grids and zones 
does need to not change during the movement of the upper 
cylinder over the lower one. By covering the cut ends of the 
cylinders with a virtual zone so that the variable solid wall 
boundary conditions become interface conditions, the up- 
per zone communicates with the lower zone only through 
the interface condition irrespective of the relative position 
of the two cylinders. As also shown in the figure, the grid 
topology of the virtual zones can be, and usually is, differ- 
ent from the topology of the volume grid in contact with 
the virtual zone. By this simple idea of converting bound- 
ary conditions into interface conditions through the use 
of virtual zones, dynamic configurations which were previ- 
ously difficult to zone and grid with patched grids become 
tractable. 

The virtual zones also allow the zonal boundaries to 
cut through the configuration surfaces, which is an im- 
portant property of control surfaces. The region of the 
configuration that intersects the zonal face is covered with 
a virtual zone to convert that region into another interface 
condition. Once a zone has been defined along with its as- 
sociated virtual zones, its definition is complete and is not 
influenced by any of its neighboring zones. In other words, 
the real zone communicates with the other zones (real or 
virtual) only through the interface conditions. Thus a par- 
ticular zone can be altered or substituted with another 
zone without any need to redefine the interface conditions 
of the other zones. For example, a zonal grid can be set 
up for a wing with control flaps with one zone for each 
(say, undeflected) flap. For the deflected flap case, only 
the flap zone needs to be replaced with a zone containing 
a deflected flap. The boundary and interface conditons of 
all other zones remain unchanged, even though they may 
now have a solid surface exposed to the flow field, e.g., the 
edge of the exposed end of the wing and flap. 

The original zoning capability of the ENSAERO code 
was extended by including the above capability of multi- 
ple interface conditions on a single block face. Since the 
code is a finite difference code the zones required an over- 
lap at the boundaries of the zones to allow for the proper 
interblock communication (i.e., interfacing). In this study 
a one- cell overlap was chosen. For this kind of zoning the 
virtual zones of zero thickness used for the finite volume 
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formulation are not appropriate. Instead, the thickness of 
the virtual zones had to be expanded to include the extent 
of the overlap of the zones and, thus, the virtual zones for 
the present formulation are now one cell thick. This formu- 
lation results in a slight mismatch of one-half cell thickness 
between the location of the actual solid wall and the vir- 
tual zone boundary. The solid wall boundary condition 
is applied at the proper location and the slight mismatch 
has no discernable influence on the overall flow field, es- 
pecially in the case where the ends of the flaps and wings 
are treated with the no-slip viscous boundary condition, 
rather than with the inviscid tangency condition. 

An example of the virtual zones required for the 
present case is shown in Figs. 4-6 for the inboard end of the 
control surface. The two virtual zones slide through each 
other with the control surface motion. It can be seen in the 
figures that different regions of the virtual zones are then 
exposed to, or in contact with, the real zones surrounding 
the wing and the flap. The virtual zones transfer the solid 
wall boundary condition to an interface condition and thus 
allows for the automatic inclusion of the variable exposure 
of the ends of the flap and wing to the flow field. The area 
where the wing and flap virtual zones overlap represents 
the unexposed portions of the flap and wing. Since the 
flow field is not updated in the virtual zones by the flow 
solver, nothing happens in the virtual zone overlap region 
nor does it influence the rest of the flow field. 


Zonal Interface Interpolation 

The original interpolation procedure used in CNSFV 
and ENSAERO is based on a global area search. Even 
though it is highly vectorized, it still requires 5000 fisec 
per target point on a Cray YMP on an interface with 2500 
points for both the target and base domains. In this pa- 
per a target point is the point to be interpolated using 
the data from the base points. The global area search is, 
thus, much too slow for a dynamic interpolation procedure 
where the interpolation coefficients have to be updated at 
every time step. By replacing the area search with a clip- 
ping search procedure based on a polygon clipping algo- 
rithm (also called window clipping, Ref. 11), the search 
time to find the interpolants is reduced by two orders of 
magnitude. This improvement brought the cpu run time 
for determining the interpolation coefficients down to the 
same level as required for the flow solver. In addition to 
the clipping search, the nested loop (or shell) search, di- 
rected hunt, and range limiting search are used to speed 
up the process and are described below. 

The clipping search is based on the well-known poly- 
gon clipping algorithm. However, instead of clipping a 
polygon to fit inside a window, only one point designated 
the target points is used. The window is formed from the 
four base points and we are trying to find the four base 
points that surround the target point. In general the four 
base points do not form a rectangular window and the 
“coarse clip” circumscribes the four base points as shown 
in Fig. 7. The clipping proceeds as follows: 

1. Circumscribe the four base points with a rough 
window. 

2. Target point to the right of left boundary? If so, 
continue with 3. If not, go to next set of base points since 
target cannot be inside this window. 


3. Target point above bottom boundary? If not, go 
to next set of base points; if so continue with 4. 

4. Target point to the left of right boundary? If not, 
go to next set of base points; if so continue with 5. 

5. Target point below top boundary? If not, go to 
next set of base points; if so target point inside the rough 
window and the “coarse clip” is successful, continue with 
6 . 

6. Transform the four base points into a unit square 
aligned with the coordinate axis using a bilinear transfor- 
mation, as sketched in Fig. 7. Locate the target point 
in the transformed space and repeat the search procedure 
2-5. If this “fine clip” is successful then have located the 
four base points that surround the target point and can 
now go to 7. 

7. Determine the bilinear interpolation coefficients 
for this target point and go to 8. 

8. Proceed to next target point and go to 1. 

The efficiency of the clipping search comes from two 
factors: the first is that non-candidate base points are re- 
jected as quickly as possible with a minimal amount of 
operations performed. Secondly, the bilinear transforma- 
tion is done only if the target point is inside the rough 
window. The bilinear transformation is relatively expen- 
sive to compute compared to setting up a rough window 
which requires only a few minmax function operations on 
the four base points, especially when the degenerate cases 
of the four base' points collapsing into triangles must, be 
considered. 

In the best-case scenario where the first guess of the 
base points is the correct one, the Cray YMP cpu time 
per target point is 35 (iscc for the above clipping search 
procedure. This includes the coarse clip, bilinear transfor- 
mation, fine clip, and determining the interpolation coef- 
ficients. The worst case is when the target points are not 
within the domain of the base points. In this case assum- 
ing grids with 2500 points for both the target and base 
grids, about 600 fisec are required. The reason that the 
times are not larger is that only the coarse clip procedure 
needs to be executed. In spite of the slow times, the worst 
case times are about 10 times faster than the vectorized 
global area search procedure. 

To obtain the times of the best case, it is clearly nec- 
essary to obtain the best guess possible for the base points 
so that only a small neighborhood needs to searched. For 
the steady-state problem a good starting guess for a new 
target point is the base points found from the previous 
target point. For unsteady cases, where the target grid is 
moving relative to the base grid, a good strategy for the 
starting base points is to use the base points found for this 
target point at the previous time step. For highly clustered 
meshes where a target point can move across several dozen 
base points during a time step, a better approach is to use 
the one for the steady state case. In this paper, however, 
we use the previous time step guess. With the starting 
guess for the base points, two procedures called the “shell 
search” and “directed hunt” are used with variable, suc- 
cess. A third method called “range limiter” is developed 
to reduce the time required for the worst case or near worst 
case situation which happens quite often for the oscillat- 
ing control surface configuration. All three methods are 
described below in further detail. 
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The shell search uses the clip search as the basic en- 
gine. The shell search refers to the sequence in which the 
four base points are tested against the clip search. The 
innermost shell depicted as ri = 1 in Fig. 8 is the starting 
guess for the base points. If the clip test is unsuccessful, 
then the search proceeds to the next shell and loops around 
the inner shell until the clip test succeeds or until all the 
8(n— 1) sets of base points are tested. If again unsuccessful, 
the search proceeds to the next shell. The number of tests 
increases as (2 n — l) 1 2 where n is the shell number. If the 
shell intersects a boundary of the base grid, then the shell 
is truncated at that boundary. The maximum number of 
shells searched is equal to the maximum dimension of the 
base grid. If the guessed starting set of base point is not 
good, then the shell search is not much better than search- 
ing the entire base grid with the clip search. However, if 
the starting guess is good and the target points have not 
moved across too many base points, then the shell search 
is successful within 2 or 3 shells, requiring at most 9 or 25 
clip searches. Tests with various cases have shown that it 
is more efficient to stop after the second shell and do the 
range limiter (explained below) before proceeding with the 
next shell. Typical test cases show that the shell hunt epu 
times vary from 35 /isec to 150 fisec per target point. 

A potentially more efficient procedure than the shell 
search is the directed hunt. In this method, the best guess 
for the four base points and the target point is transformed 
into the unit square window. Let the transformed coordi- 
nates of the target point be ( k tl l t ) as shown in Fig. 9. 
If the coordinates are positive and less than one, the four 
base points have been found; if not, the next guess for the 
reference base point coordinates is given by 


k n = k Q -j- int(k t — 1) 

Jn = J© + int(l t - 1) 

where k Q Jd are the indices of the reference point of the 
four base points, as shown in Fig. 9. Typically, the di- 
rected hunt should find the base points with only 2 to 3 
steps. However the bilinear transformation often breaks 
down because of degenerate base points, or the direction 
takes the search to the outside of the base domain bound- 
aries (e.g., in crossing the wake boundary of a C-mesh). 
In the latter case one has to step along the boundary until 
the directed hunt can be continued. For these reasons the 
directed hunt has not been reliable and efficient enough 
for the oscillating control surface problem. Since the di- 
rected hunt has the potential of being 2 to 3 times faster 
than the shell search, especially if the first guess is poor, 
work is continuing on developing a reliable directed hunt 
procedure. 

The cost of computing the searching procedure in- 
creases as the number of base points increases. What is 
needed is an efficient procedure to reduce the number of 
base points over which to search, or to eliminate them en- 
tirely if the target point is not within the domain of the 
i base points as in the worst-case scenario mentioned above. 
The “range limiter” was developed for this purpose. As 
with the shell search the basic engine is again the clip test. 
The range limiter can be described with the help of Fig. 
10 as follows: 

1. Circumscribe the domain of the base points with 

the window shown in Fig. 10a. 


2. Apply the clip test to the target point; if target is 
inside window, continue with 3; if not. target point cannot 
be interpolated with the given base points. Go on to next 
target point. 

3. Subdivide the base domain in half and choose the 
subdivision which produces the smallest overlap area as 
shown in Fig. 10b and 10c. If overlap areas axe equal, 
choose the partition which operates on the larger of the 
two dimensions; if the dimensions are the same, use the 
first of the two subdivisions. 

4. Apply the clip test to the target point with both 
of the subdomains. If the target point is in one but not 
the other of the subdomains, then define the domain of 
base points to be the domain in which the target point lies 
and continue with 1. If the target point lies in both of 
the subdomains (i.e., it lies in the overlap area) then cut 
off the two subdomains at the ends opposite the overlap 
region and redefine the domain to be the central portion 
of the base point domain and continue with 1. The clip 
test is applied to the cutoff portions to make sure that the 
trimmed portions do not contain the target point. 

5. The subdivision along any one direction stops if 
the dimension in that direction is less than 5, since at that 
point the shell search is more efficient. 

Although the range limiter appears cumbersome, it 
is quite useful and efficient especially if the domain of base 
points is larger than the domain of target points, and it 
very quickly excises the portions of the base domain which 
are not close to the target point. Once the base domain 
has been trimmed, the shell search or directed hunt is used 
for the local search. Another way of trimming the base 
domain is to place search limits on the base domain beyond 
which the search will not be conducted. This last method 
is effective only if the search limits are known a priori, as 
for the present case of the oscillating control surface. 

It is of interest to compare the present method with 
the so-called Domain Connectivity Function method of 
Ref. 12. The present method is just slightly faster than 
the dynamic mode of the DCF method. Compared to the 
static mode of the DCF method, where the inverse map- 
ping needs to be determined, the present method is sub- 
stantially faster (approx. 10 times). For the oscillating 
control surface problem, the DCF method would have to 
be run in the static mode and thus is too slow to be prat- 
ical. 

Results 

The test case considered in the present study is a 
clipped delta wing with an oscillating trailing edge control 
surface [13]. The wing planform is shown in Fig. 1. The 
wing has a leading edge sweep angle of 50.4 deg and a 
6% thick circular arc airfoil section. At = 0.9 and 
a = 3 deg, both a leading edge vortex and a shock wave are 
present on the upper surface of the wing. The C-H grids of 
the three real zones consist of 151 x 13 x 34, 151 x 15 x 34, 
and 151 x 20 x 34 points from inboard to outboard as shown 
in Fig. 2. Since the experiment was conducted using a 
Freon test medium, the ratio of specific heats, 7, is set 
to 1.135 in the present computations. As stated before, 
the modified Baldwin-Lomax model is used to account for 
the leading edge and control surface vortices. However, 
the boundaries normal to the edges of the control surface 


H 



and the wing cuts are treated with the solid wall tangency 
condition. Steady state and rigid pitching calculations of 
this wing were reported in Ref. 8. 

Figure 11 shows the unsteady pressures coefficients 
with the control surface oscillating at a frequency of 8 Hz 
and an amplitude of 6.65 deg at Moo = 0.9, a = 3 deg 
and Re c = 17 x 10 6 based on the root chord. The results 
are shown as the amplitude and phase angle of the upper 
surface pressure coefficients at the three span stations as 
indicated on the figure. In general, the agreement with the 
experimental results is good. Because the accuracy of ex- 
periment has its own limitations, the virtual zone results 
are also compared with the single grid results of Ref. 5. 
There is quite a discrepancy between the virtual zone re- 
sults and the single grid results (the 151 x 44 x 34 grid), 
especially in the amplitudes at the center of the control 
surface. This discrepancy is most likely due to the gap 
that was introduced between the flap and wing in the sin- 
gle grid case to accommodate the shearing grid. If the 
gap is reduced by increasing the spanwise resolution of the 
wing and control surface, the computational results of the 
refined single grid case (151 X 87 x 34) approach those of 
the virtual zone. This particular example demonstrates 
the importance of simulating the geometry of the control 
surface/ wing configuration accurately. 

Figure 12 shows the upper surface pressures as well 
as the instantaneous streamline traces emanating from the 
leading edge of the wing (black traces), the lower edges of 
the wing at the control surface cut (blue traces), and the 
upper edges of the control flaps (red traces). The inter- 
action between the leading edge vortex and the outboard 
edge of the control surface vortices during the flap motion 
cycle is well demonstrated. Although not apparent in this 
figure, the computations also show that the flow on the 
leading side (defined as upper surface during the upstroke 
motion and lower surface during the downstroke motion) 
of the inboard section of the control surface is separated. 

Figure 13 shows the velocity vectors at the inboard 
edge of the control surface as the flap is in the upstroke 
mode at the instant of zero deflection. The flow field is 
shown in much more detail at the trailing edge of the con- 
trol surface. The detailed views show the velocities of 
the trailing edge of the control surface. As can be seen, 
these velocities are small, even when compared to the flow 
field in the inner regions of the boundary layer. Figure 14 
shows the surface streamlines on the inboard edge and up- 
per surface of the control surface. The particle traces are 
computed at each instant of time by freezing the flow field 
during the tracing procedure. The temporal changes in 
the surface particle traces as the control surface oscillation 
cycle is executed are quite apparent, both on the upper 
surface and the edge of the control surface as it is exposed 
and covered by the wing edge. The first view of Figure 
14 shows that the flow is separated near the inboard trail- 
ing edge corner of the control surface. Preliminary studies 
with a viscous (albeit laminar) treatment of the control 
surface edges and wing cuts indicate that the separated 
region on the upper surface of the flap increases. 

Conclusions 

An unsteady interface algorithm based on the idea 
of virtual zones has been developed for a finite difference 
code, ENSAERO. The new “virtual” zoning technique sim- 
plifies the zoning of complex geometries, such as control 


surfaces, and makes possible the use of standard multizonal 
codes for complex configurations. A fast search routine 
based on a window clipping algorithm has also been devel- 
oped and is fast enough for the interpolation coefficients 
to be recomputed at every time step. For the example pre- 
sented in this study, the computational effort required for 
the interpolation coefficients is less than that for the flow 
solver. 

Both of the above developments have made practical 
a complete unsteady Navier- Stokes simulation of a forced 
oscillating control surface on a clipped delta wing in the 
transonic flow regime. The method has been validated 
against experimental data as well as numerical simulations 
based on a shearing single-zone grid. The numerical result 
confirms that the accurate representation of geometry by 
“virtual zones” is superior to the single zone computations 
of the same grid size. 
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Fig. 2. Surface grids and zonal boundaries of 
wing/flap configuration. 


Fig. 1. Planform of clipped delta wing with trailing 
edge flap. 



(i) Conventional Grid OK 



( ii ) Conventional Grid ?? 



( iii ) Conventional Grid OK 


Fig. 3a. Generic wing/flap configuration with con- 
ventional zones. The top cylinder is sliding over the face 
of the bottom cylinder. 



Fig. 3l>. Generic wing/flap configuration demon- 
strating the virtual zone concept. The top cylinder is slid- 
ing over rhe face of the bottom cylinder. 



Fig. 4. Perspective view of trailing edge flap and 

wing. 



Fig. 5. Perspective view of trailing edge flap and 
wing with the inboard wing and flap virtual zones. 



Fig. 6. Perspective phantom view of trailing edge 
flap and wing with the inboard wing and flap virtual zones. 
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Fig. 7. Clip search procedure; Coarse clip, bilinear 
transformation and fine clip test. 


Fig. 9. The directed hunt procedure. 
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(a) 


Fig. 10. The range limiter procedure. 


Fig. 8. The sequence of steps required for the shell 
search procedure. 
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Fig. 10. The range limiter procedure. 
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Fig. 11. Comparison of unsteady pressure coeffi- 
cients between experiment (Ref. 13) and computations 
using virtual zones and single grids. 









Fig. 12. Upper surface pressures and instantaneous 
streamlines emanating from the leading edge and the edges 
of the control surface. 
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Fig. 13. Velocity vector plots near the inboard edge 
of the control surface at the instant of zero deflection and 
maximum flap velocity. 
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Fig. 14. Surface particles traces on the upper surface 
and inboard edge of the control surface at four instances 
of the control surface cycle. 
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